pardeu Subroutine

      pure subroutine pardeu(tx,nx,ty,ny,c,kx,ky,nux,nuy,x,y,z,m,wrk,lwrk,iwrk,kwrk,ier)

      !  subroutine pardeu evaluates on a set of points (x(i),y(i)),i=1,...,m the partial derivative
      !  ( order nux,nuy) of a bivariate spline s(x,y) of degrees kx and ky, given in the b-spline
      !  representation.
      !
      !  calling sequence:
      !     call pardeu(tx,nx,ty,ny,c,kx,ky,nux,nuy,x,y,z,m,wrk,lwrk,iwrk,kwrk,ier)
      !
      !  input parameters:
      !   tx    : real array, length nx, which contains the position of the knots in the x-direction.
      !   nx    : integer, giving the total number of knots in the x-direction
      !   ty    : real array, length ny, which contains the position of the knots in the y-direction.
      !   ny    : integer, giving the total number of knots in the y-direction
      !   c     : real array, length (nx-kx-1)*(ny-ky-1), which contains the b-spline coefficients.
      !   kx,ky : integer values, giving the degrees of the spline.
      !   nux/y : integer values, specifying the order of the partial derivative. 0<=nux<kx, 0<=nuy<ky.
      !   kx,ky : integer values, giving the degrees of the spline.
      !   x     : real array of dimension (mx).
      !   y     : real array of dimension (my).
      !   m     : on entry m must specify the number points. m >= 1.
      !   wrk   : real array of dimension lwrk. used as workspace.
      !   lwrk  : integer, specifying the dimension of wrk.
      !           lwrk >= mx*(kx+1-nux)+my*(ky+1-nuy)+(nx-kx-1)*(ny-ky-1)
      !   iwrk  : integer array of dimension kwrk. used as workspace.
      !   kwrk  : integer, specifying the dimension of iwrk. kwrk >= mx+my.
      !
      !  output parameters:
      !   z     : real array of dimension (m).
      !           on successful exit z(i) contains the value of the specified partial derivative of s(x,y)
      !           at the point (x(i),y(i)),i=1,...,m.
      !   ier   : integer error flag
      !
      !  restrictions:
      !   lwrk>=m*(kx+1-nux)+m*(ky+1-nuy)+(nx-kx-1)*(ny-ky-1),
      !
      !  other subroutines required:
      !    fpbisp,fpbspl
      !
      !  references :
      !    de boor c : on calculating with b-splines, j. approximation theory 6 (1972) 50-62.
      !   dierckx p. : curve and surface fitting with splines, oxford university press, 1993.
      !
      !  author :
      !    p.dierckx
      !    dept. computer science, k.u.leuven
      !    celestijnenlaan 200a, b-3001 heverlee, belgium.
      !    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
      !
      !
      !  ..scalar arguments..
      integer    , intent(in)    :: nx,ny,kx,ky,m,lwrk,kwrk,nux,nuy
      integer    , intent(out)   :: ier
      !  ..array arguments..
      integer    , intent(inout) :: iwrk(kwrk)
      real(RKIND), intent(in)    :: tx(nx),ty(ny),c((nx-kx-1)*(ny-ky-1)),x(m),y(m)
      real(RKIND), intent(out)   :: z(m)
      real(RKIND), intent(inout) :: wrk(lwrk)

      !  ..local scalars..
      integer :: i,iwx,iwy,j,kkx,kky,kx1,ky1,lx,ly,lwest,l1,l2,mm,m0,m1,nc,nkx1,nky1,nxx,nyy
      real(RKIND) :: ak,fac

      !  ..
      !  before starting computations a data check is made. if the input data are invalid control is
      !  immediately repassed to the calling program.
      ier   = FITPACK_INPUT_ERROR
      kx1   = kx+1
      ky1   = ky+1
      nkx1  = nx-kx1
      nky1  = ny-ky1
      nc    = nkx1*nky1
      lwest = nc +(kx1-nux)*m+(ky1-nuy)*m
      if (nux<0 .or. nux>=kx) return
      if (nuy<0 .or. nuy>=ky) return
      if (lwrk<lwest)         return
      if (kwrk<(m+m))         return
      if (m<1)                return

      ier = FITPACK_OK
      nxx = nkx1
      nyy = nky1
      kkx = kx
      kky = ky

      !  the partial derivative of order (nux,nuy) of a bivariate spline of degrees kx,ky is a bivariate
      !  spline of degrees kx-nux,ky-nuy. we calculate the b-spline coefficients of this spline
      wrk(:nc) = c(:nc)

      if (nux>0) then
          lx = 1
          x_deriv_order: do j=1,nux
             ak = kkx
             nxx = nxx-1
             l1 = lx
             m0 = 1
             do i=1,nxx
                l1 = l1+1
                l2 = l1+kkx
                fac = tx(l2)-tx(l1)
                if (fac>zero) THEN
                   do mm=1,nyy
                      m1 = m0+nyy
                      wrk(m0) = (wrk(m1)-wrk(m0))*ak/fac
                      m0  = m0+1
                   end do
                endif
             end do
             lx = lx+1
             kkx = kkx-1
         end do x_deriv_order
      endif

      if (nuy>0) then
          ly = 1
          y_deriv_order: do j=1,nuy
             ak = kky
             nyy = nyy-1
             l1 = ly
             do i=1,nyy
                l1 = l1+1
                l2 = l1+kky
                fac = ty(l2)-ty(l1)
                if (fac>zero) then
                   m0 = i
                   do mm=1,nxx
                      m1 = m0+1
                      wrk(m0) = (wrk(m1)-wrk(m0))*ak/fac
                      m0  = m0+nky1
                   end do
                endif
             end do
             ly = ly+1
             kky = kky-1
          end do y_deriv_order
          m0 = nyy
          m1 = nky1
          do mm=2,nxx
            do i=1,nyy
              m0 = m0+1
              m1 = m1+1
              wrk(m0) = wrk(m1)
            end do
            m1 = m1+nuy
          end do
      endif

      !  we partition the working space and evaluate the partial derivative
      iwx = 1+nxx*nyy
      iwy = iwx+m*(kx1-nux)

      do i=1,m
         call fpbisp(tx(nux+1),nx-2*nux,ty(nuy+1),ny-2*nuy,wrk,kkx,kky, &
                     x(i),1,y(i),1,z(i),wrk(iwx),wrk(iwy),iwrk(1),iwrk(2))
      end do
      return
      end subroutine pardeu